pacman::p_load(sf, raster, spatstat, tmap, tidyverse,sparr,spNetwork,dplyr,animation,stringr)Take-home_Ex03
Take home exe 03
Data
Philippines data from Armed Conflict Location & Event Data Project (ACLED).
Philippines boundary data from HDX
Importing data and package
package
- spNetwork, which provides functions to perform Spatial Point Patterns Analysis such as kernel density estimation (KDE) and K-function on network. It also can be used to build spatial matrices (‘listw’ objects like in ‘spdep’ package) to conduct any kind of traditional spatial analysis with spatial weights based on reticular distances.
- sf package provides functions to manage, processing, and manipulate Simple Features, a formal geospatial data standard that specifies a storage and access model of spatial geometries such as points, lines, and polygons.
- tmap which provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API.
- spatstat, which has a wide range of useful functions for point pattern analysis. In this hands-on exercise, it will be used to perform 1st- and 2nd-order spatial point patterns analysis and derive kernel density estimation (KDE) layer.
- raster which reads, writes, manipulates, analyses and model of gridded spatial data (i.e. raster). In this hands-on exercise, it will be used to convert image output generate by spatstat into raster format.
- maptools which provides a set of tools for manipulating geographic data. In this hands-on exercise, we mainly use it to convert Spatial objects into ppp format of spatstat.
data
drug_case <- read_csv("data/2016-01-01-2024-06-30-Philippines.csv")%>%
st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
st_transform(crs = 32651) %>%
mutate(event_date = dmy(event_date)) %>%
mutate(event_month = year*100 + month(event_date)) %>%
mutate(event_quarter = year*10 + quarter(event_date))%>%
mutate(quarter = as.numeric(str_sub(event_quarter, 5, 5)))Rows: 6921 Columns: 31
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (23): event_id_cnty, event_date, disorder_type, event_type, sub_event_ty...
dbl (8): year, time_precision, iso, latitude, longitude, geo_precision, fat...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.import the boundary data
ph_sf = st_read(dsn = "data/phl_adm_psa_namria_20231106_shp",layer = "phl_admbnda_adm2_psa_namria_20231106")Reading layer `phl_admbnda_adm2_psa_namria_20231106' from data source
`/Users/liangyuhang/Downloads/Maaaaaaaaaark/IS415_g/Take-home_Ex/Take-home_Ex03/data/phl_adm_psa_namria_20231106_shp'
using driver `ESRI Shapefile'
Simple feature collection with 88 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 114.2779 ymin: 4.587294 xmax: 126.605 ymax: 21.12189
Geodetic CRS: WGS 84chack the crs
st_crs(ph_sf)Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]transform the crs
ph_sf = st_transform(ph_sf,crs = 32651)
st_crs(ph_sf)Coordinate Reference System:
User input: EPSG:32651
wkt:
PROJCRS["WGS 84 / UTM zone 51N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 51N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",123,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["Between 120°E and 126°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Japan. North Korea. Philippines. Russian Federation. South Korea. Taiwan."],
BBOX[0,120,84,126]],
ID["EPSG",32651]]ph_sf_owin <-as.owin(ph_sf)show the point
tmap_mode('view')tmap mode set to interactive viewingtm_shape(drug_case)+tm_dots()